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Deterministic treatment of model error in 
geophysical data assimilation 


Alberto Carrassi and Stephane Vannitsem 


Abstract This chapter describes a novel approach for the treatment of model error in 
geophysical data assimilation. In this method, model error is treated as a determin¬ 
istic process fully correlated in time. This allows for the derivation of the evolution 
equations for the relevant moments of the model error statistics required in data as¬ 
similation procedures, along with an approximation suitable for application to large 
numerical models typical of environmental science. In this contribution we hrst de¬ 
rive the equations for the model error dynamics in the general case, and then for the 
particular situation of parametric error. We show how this deterministic description 
of the model error can be incorporated in sequential and variational data assimi¬ 
lation procedures. A numerical comparison with standard methods is given using 
low-order dynamical systems, prototypes of atmospheric circulation, and a realis¬ 
tic soil model. The deterministic approach proves to be very competitive with only 
minor additional computational cost. Most importantly, it offers a new way to ad¬ 
dress the problem of accounting for model error in data assimilation that can easily 
be implemented in systems of increasing complexity and in the context of modern 
ensemble-based procedures. 


1 Introduction 


The prediction problem in geophysical fluid dynamics typically relies on two com¬ 
plementary elements: the model and the data. The mathematical model, and its dis¬ 
cretized version, embodies our knowledge about the laws governing the system evo¬ 
lution, while the data are samples of the system’s state. They give complementary 
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information about the same object. The sequence of operations that merges model 
and data to obtain a possibly improved estimate of the flows state is usually known, 
in environmental science, as data assimilation Hi ED. The physical and dynamical 
complexity of geophysical systems makes the data assimilation problem particularly 
involved. 

The different information entering the data assimilation procedure, usually the 
model, the data and a background field representing the state estimate prior to the 
assimilation of new observations, are weighted according to their respective accu¬ 
racy. Data assimilation in geophysics, particularly in numerical weather prediction 
(NWP) has experienced a long and fruitful stream of research in recent decades 
which has led to a number of advanced methods able to take full advantage of the 
increasing amount of available observations and to efficiently track and reduce the 
dynamical instabilities m. As a result the overall accuracy of the Earths system 
estimate and prediction, particularly the atmosphere, has improved dramatically. 

Despite this trend of improvement, the treatment of model error in data assimila¬ 
tion procedures is still, in most instances, done following simple assumptions such 
as the absence of time correlation ifT^ . The lack of attention on model error is in part 
justified by the fact that on the time scale of NWP, where most of the geophysical 
data assimilation advancements have been originally concentrated, its influence is 
reasonably considered small as compared to the initial condition error that grows in 
view of the chaotic nature of the dynamics. Nevertheless, the improvement in data 
assimilation techniques and observational networks on the one hand, and the recent 
growth of interest in seasonal-to-decadal prediction on the other lfT3ll47l . has placed 
model error, and its treatment in data assimilation, as a main concern and a key pri¬ 
ority. A number of studies reflecting this concern have appeared, in the context of 
sequential and variational schemes ifT^l42l l43ll23l . 

Two main obstacles toward the development of techniques taking into account 
model error sources are the huge size of the geophysical models and the wide range 
of possible model error sources. The former problem implies the need to estimate 
large error covariance matrices on the basis of the limited number of available ob¬ 
servations. The second important issue is related to the multiple sources of modeling 
error, such as incorrect parametrisation, numerical discretization, and the lack of de¬ 
scription of some relevant scale of motion. This latter problem has until recently lim¬ 
ited the development of a general formulation for the model error dynamics. Model 
error is commonly modeled as an additive, stationary, zero-centered, Gaussian white 
noise process. This choice could be legitimate by the multitude of unknown error 
sources and the central limit theorem. However, despite this simplification, the size 
of geoscientific models still makes detailed estimation of the stochastic model error 
covariance impractical. 

In the present contribution we describe an alternative approach in which the evo¬ 
lution of the model error is described based on a deterministic short-time approxi¬ 
mation. The approximation is suitable for realistic applications and is used to esti¬ 
mate the model error contribution in the state estimate. The method is based on the 
theory of deterministic dynamics of the model error that was introduced recently by 
EalsiEa- Using this approach it is possible to derive evolution equations for the 


Deterministic treatment of model error in geophysical data assimilation 


3 


moments of the model error statistics required in data assimilation procedures, and 
has been applied in the context of both sequential and variational data assimilation 
schemes, and for errors originated from uncertain parameters and from unresolved 
scales. 

We give here a review of the recent developments of the deterministic treatment 
of model error in data assimilation. To this end, we start by first formalizing the 
deterministic model error dynamics in Sect. 2. We show how general equations for 
the mean and covariance error can be obtained and discuss the parametric error as a 
special case. In Sections 3 and 4 the incorporation of the short-time model error evo¬ 
lution laws is described in the context of the Extended Kalman filter and variational 
scheme respectively. These two types of assimilation procedures are significantly 
different and are summarized in the respective Sections along with the discussion on 
the consequences of the implementation of the model error treatment. We provide 
some numerical illustrations of the proposed approach together with comparisons 
with other methods, for two prototypical low order chaotic systems widely used in 
theoretical studies in geosciences Il28ll29l and a quasi-operational soil model ll^ . 

New potential applications of the use of the deterministic model error treatment 
are currently under way and are summarized, along with a synopsis of the method, 
in the final discussion Section 5. These include soil data assimilation with the use 
of new observations and ensemble based procedures na. 


2 Formulation 


Let the model at our disposal be represented as: 


dx{t) 

dt 


f(x,A), 


( 1 ) 


where f is typically a nonlinear function, defined in and A is a P-dimensional 
vector of parameters. 

Model ([T]i is used to describe the evolution of a (unknown) true dynamics, i.e. 
nature, whose evolution is assumed to be given by the following coupled equations: 


dx{t) 

dt 


f(x,y,A) 


dyjt) 

dt 


= g(x,y,A) 


( 2 ) 


where x is a vector in R^, and y is defined in R^ and may represent scales that are 
present in the real world, but are neglected in model ([T]i; the unknown parameters 
A have dimension P. The true state is thus a vector of dimension N + L. The model 
state vector x and the variable x of the true dynamics span the same phase space 
although, given the difference in the functions f and f, they do not have the same 
attractor in general. The function f can have an explicit dependence on time but it is 
dropped here to simplify the notation. 
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When using model ([TJ to describe the evolution of x, estimation error can arise 
from the uncertainty in the initial conditions at the resolved scale (x(fo) 7^ x(fo)) 
and from the approximate description of the nature afforded by ([T]i which is referred 
as model error. A number of different sources of model errors are present in en¬ 
vironmental modeling. Typical examples are those arising from the inadequate de¬ 
scription of some physical processes, numerical discretization and/or the presence 
of scales in the actual dynamics that are unresolved by the model. The latter are 
typically parametrised in terms of the resolved variables (for instance the Reynolds 
stress of the turbulent transport). 


2.1 General description of model error dynamics 

Following the approach outlined in 041 . we derive the evolution equations of the 
dominant moments, mean and covariance, of the estimation error 5x = x — x in the 
resolved scale {i.e. in R*^). The formal solutions of ([T]i and (|2|i read respectively: 



(3) 



(4) 


where xq = x(fo), and xq = x(fo). By taking the difference between © and (HI, and 
averaging over an ensemble of perturbations around a reference state, we get the 
formal solution for the mean error, the bias: 



(5) 


with 5xo = xo — Xq. Two types of averaging could be performed, one over a set of 
initial conditions sampled on the attractor of the system, and/or a set of perturbations 
around one specific initial state selected on the system’s attractor. In data assimila¬ 
tion, the second is more relevant since one is interested in the local evaluation of the 
uncertainty. However, in many situations the first one is used to get statistical infor¬ 
mation on covariances quantities, as will be illustrated in this Chapter. For clarity, 
we will refer to < . > as the local averaging, and to << . >> for an averaging over 
a set of initial conditions sampled over the attractor of the system. In this section, 
we will only use < . > for clarity, but it also extends to the other averaging. We will 
use the other notation << . >> when necessary. 

In the hypothesis that the initial condition is unbiased, < 5xo >= 0, Eq. Q gives 
the evolution equation of the bias due to the model error, usually refers to as drift 
in climate prediction context. The important factor driving the drift is the difference 
between the true and modeled tendency fields, < f(x(T),A) -f(x(T),y(T),A) >. 
Expanding Q in Taylor series around fo = 0 up to the first non-trivial order, and 
using unbiased initial conditions, it reads: 
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b'"(f) =< 5x(f) >Ri<f(x(T),A)-f(x(T),y(T),A) > t (6) 

Equation (|6]l gives the evolution of the bias, b™, the drift, in the short-time approx¬ 
imation and the subscript m stands for model error-related bias. It is important to 
remark that in the case of stochastic model error treatment, and in the hypothesis of 
unbiased initial condition error, b™ = 0. 

Similarly, by taking the expectation of the external product of the error anomalies 
5x by themselves, we have; 

P(f) =< {5x(f)}{5x(f)}^ >=< {5xo}{5xo}^ > + 

< {5xo}{ j|'r/T[f(x(T),A)-f(x(T),y(T),A)]}^ > + 

< {jJ'r/T[f(x(T),A) -f(x(T),y(T),A)]}{5xo}^ > + 

[ dz [ dz < {f(x(T),A)-f(x(T),y(T),A)}{f(x(T'),A)-f(x(T'),y(T'),A)}^ > 
Jo Jo 

(7) 

Equation (|7]) describes the time evolution of the estimation error covariance in the 
resolved scale. The first term, that does not depend on time, represents the covari¬ 
ance of the initial error. The two following terms account for the correlation between 
the error in the initial condition and the model error, while the last term combines 
the effect of both errors on the evolution of the estimation error covariance. 

Let us focus on the last term of Eq. (I7]i denoted as. 


P(r) = 


JO JO 


dz < {f(x(T),A)-f(x(T),y(T),A)}{f(x(T'),A)-f(x(T'),y(T'),A)}^ 


( 8 ) 

The amplitude and structure of this covariance depends on the dynamical proper¬ 
ties of the difference of the nature and model tendency fields. Assuming that these 
differences are correlated in time, we can expand ® in a time series up to the first 
nontrivial order around the arbitrary initial time to = 0, and gets: 


P'"(f) Ri< {f(xo,A) -f(xo,yo,A)}{f(xo,A) -f(xo,yo,A)}^ > = Qf^ (9) 


where Q is the model error covariance matrix at initial time. Note again that, if the 
terms f — f are represented as white-noise process, the short-time evolution of P(f) 
is bound to be linear instead of quadratic. This distinctive feature is relevant in data 
assimilation applications where model error is often assumed to be uncorrelated in 
time, a choice allowing for a reduction of the computational cost associated with 
certain types of algorithms iiiia. 
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2,2 Model error due to parameter uncertainties 


We assume for simplicity that the model resolves all scales present in the reference 
system. Under the aforementioned hypothesis that the model and the true trajectories 
span the same phase space, nature dynamics, can be rewritten as: 

=f(x,A) + eh(x, 7 ) ( 10 ) 

The function h, which has the same order of magnitude of f and is scaled by the 
dimensionless parameter e, accounts for all other extra terms not included in the 
model and depends on the resolved variable x and on a set of additional parameters 
7 . In a more formal description, this h would correspond to a function relating the 
variables x and y under an adiabatic elimination ll34l . We are interested here in a 
situation in which the main component of the nature dynamics is well captured by 
the model so that e << 1, and the extra terms described by h are neglected. We 
concentrate in a situation in which model error is due only to uncertainties in the 
specification of the parameters appearing in the evolution law f. This formulation 
accounts, for instance, for errors in the description of some physical processes (dis¬ 
sipation, external forcing, etc.) represented by the parameters. 

An equation for the evolution of the state estimation error 5x can be obtained 
by taking the difference between the first rhs term in (flOl i and ([T]i. The evolution 
of 5x depends on the error estimate at the initial time f = fo (initial condition error 
5x(fo) = 5xo) and on the model error. If 5x is ’’small”, the linearized dynamics 
provides a reliable approximation of the actual error evolution. The linearization is 
made along a model trajectory, solution of ([ 1 ]), by expanding, to first order in 5x and 
5A = A — A, the difference between Eqs. (fTOl i and ([T]i: 


d5x 

dt 


-I 

5x' 


.5x- 


-I 

dX' 


i5A 


( 11 ) 


The first partial derivative on the rhs of (fTTT i is the Jacobian of the model dynamics 
evaluated along its trajectory. The second term, which corresponds to the model 
error, will be denoted 5jU hereafter to simplify the notation; 5p. = 

The solution of (fTTIi . with initial condition 5xo at f = fo, reads: 


5x{t) = M,,,^5xo + 


f 

Jt(\ 


dT:Mt,T:5iJ.{z) 


= 5x'"(f) + 5x'''(f) ( 12 ) 

with M, being the fundamental matrix (the propagator) relative to the linearized 
dynamics along the trajectory between to and f. We point out that 5p and Mr.^ in 
(fT^ depend on T (the integration variable) through the state variable x. Equation 
([Ell states that, in the linear approximation, the error in the state estimate is given 
by the sum of two terms, the evolution of initial condition error, 5x“^, and the model 
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error, Sx™. The presence of the fundamental matrix M in the expression for 5x"’ 
suggests that the instabilities of the flow plays a role in the dynamics of model error. 

Let us now apply the expectation operator to (fT2li defined locally around the 
reference trajectory, by sampling over an ensemble of initial conditions and model 
errors, and the equation for the mean estimation error along a reference trajectory 
reads: ^ 

< dx{t) >=M,tQ < dxo> + diM, t: < dn {T) > 

JtQ 

=< 5x'" > + < ax'" > (13) 

In a perfect model scenario an unbiased state estimate at time to (< ^xq >= 0) 
will evolve, under the linearized dynamics, into an unbiased estimate at time t. In 
the presence of model error and, depending on its properties, an initially unbiased 
estimate can evolve into a biased one with < 5n{t) > being the key factor. 

The dynamics of the state estimation error covariance matrix can be obtained by 
taking the expectation of the outer product of 5x{t) with itself. Assuming that the 
estimation error bias is known and removed from the background error, we get: 

P(f) =< 5x{t)5x{t)^ > 

= P'"(f) +P'’'(f) +P“"(0 + (P“'''')^(0 (14) 

where: 

P'‘(f) = M,,,o < 5xo5xo^ > (15) 

P'"(f)= (16) 

Jto Jto 

P™'-'-(f) = < (5xo) (^J‘ dTM,^,5^liT)^ > (17) 

The four terms of the r.h.s. of (fT4li give the contribution to the estimation error co- 
variance at time t due to the initial condition, model error and their cross correlation, 
respectively. These integral equations are of little practical use for any realistic non¬ 
linear systems, let alone the big models used in environmental prediction. A suitable 
expression can be obtained by considering their short-time approximations through 
a Taylor expansion around Iq. We proceed by expanding (fT2l i in Taylor series, up to 
the first non trivial order, only for the model error term 5x™ while keeping the initial 
condition term, Sx''^, unchanged. In this case, the model error 5x'" evolves linearly 
with time according to: 

5x"'Ri 5jUo(f-fo) (18) 

where 5jU(fo) = 

By adding the initial condition error term, 8x“^, we get a short time approxima¬ 
tion of (fT2l) : 

5x(f) R! M,,;q5xo- f 5jUo(f - fo) 

For the mean error we get: 


(19) 
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b'"(f) «< 5x(f) >« < 5x0 > + < 5/io >{t- to) (20) 

Therefore, as long as < 5/Xo > is different from zero, the bias due to parametric er¬ 
ror evolves linearly for short-time, otherwise the evolution is conditioned by higher 
orders of the Taylor expansion. Note that the two terms in the short time error evo¬ 
lution (fT^ and (l20l i. are not on equal footing since, in contrast to the model error 
term, which has been expanded up to the first nontrivial order in time, the initial 
condition error evolution contains all the orders of times The point is 

that, as explained below, we intend to use these equations to model the error evolu¬ 
tion in conjunction with the technique of data assimilation for which the full matrix 
M, or an amended ensemble based approximation, is already available. 

Taking the expectation value of the external product of (fT^ by itself and averag¬ 
ing, we get: 

P(f) « < 5xo5xo^ > + 

-f [< 5/io5xo^ < dxod^lf/' >](f-fo)+ < S^lo^n/ > {t-tof 

( 21 ) 

Equation (l2Tli is the short time evolution equation, in this linearized setting, for 
the error covariance matrix in the presence of both initial condition and parametric 
model errors. 


3 Deterministic model error treatment in the extended Kalman 
filter 

We describe here two formulations of the extended Kalman filter (EKE) incorporat¬ 
ing a model error treatment. The Short-Time-Extended-Kalman-Eilter, ST-EKE ID 
accounts for model error through an estimate of its contribution to the assumed 
forecast error statistics. In the second formulation, the Short-Time-Augmented- 
Extended-Kalman-Eilter, ST-AEKE ||6|, the state estimation in the EKE is accompa¬ 
nied with the estimation of the uncertain parameters. This is done in the context of 
a general framework known as state augmentation 03 . In both cases model error is 
treated as a deterministic process implying that the dynamical laws described in the 
previous section are incorporated, at different stages, in the filter formulations. 

The EKE extends, to nonlinear dynamics, the classical Kalman filter (KE) for 
linear dynamics ll20l . The algorithm is sequential in the sense that a prediction of 
the system’s state is updated at discrete times, when observations are present. The 
state update, the analysis, is then taken as the initial condition for the subsequent 
prediction up to the next observation time. The EKE, as well as the standard KE for 
linear dynamics, is derived in the hypothesis of Gaussian errors whose distributions 
can thus be fully described using only the first two moments, the mean and the 
covariance. Although this can represent a very crude approximation, especially for 
nonlinear systems, it allows for a dramatic reduction of the cost and difficulties 
involved in the time propagation of the full error distribution. 
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The model equations can conveniently be written in terms of a discrete mapping 
from time 4 to tt+i- 

x(^i=^x^ (22) 

where x-^ and x“ are the forecast and analysis states respectively and ^ is the non¬ 
linear model forward operator (the resolvent of ([T]i). 

Let us assume that a set of M noisy observations of the true system (l2]i, stored 
as the components of an M-dimensional observation vector y°, is available at the 
regularly spaced discrete times ti^ = tQ + kx,k=\,2..., with T being the assimilation 
interval, so that; 

y^' = ,^(x,) + e^ (23) 

where e” is the observation error, assumed here to be Gaussian with known covari¬ 
ance matrix R and uncorrelated in time. M’ is the (possibly nonlinear) observation 
operator which maps from model to observation space {i.e. from model to observed 
variables) and may involve spatial interpolations as well as transformations based 
on physical laws for indirect measurements ifTSll . 

For the EKF, as well as for most least-square based assimilation schemes, the 
analysis state update equation at an arbitrary analysis time 4 , reads Gl: 

x" = [I - KH] Ttf + Ky" (24) 

where the time indexes are dropped to simplify the notation. The analysis error 
covariance, P", is updated through: 

P" = [I-KH]P^ (25) 

The I xM gain matrix K is given by: 

K = P^ [HP^H^ + R] ^ ’ (26) 

where P^ is the I x I forecast error covariance matrix and H the linearized observa¬ 
tion operator {aMxI real matrix). The analysis update is thus based on two comple¬ 
mentary sources of information, the observations, y", and the forecast x^. The errors 
associated to each of them are assumed to be uncorrelated and fully described by 
the covariance matrices R and P-^ respectively. 

In the EKF, the forecast error covariance matrix, P^, is obtained by linearizing 
the model around its trajectory between two successive analysis times 4 and = 
fjt -f T. In the standard formulation of the EKF model error is assumed to be a random 
uncorrelated noise whose effect is modeled by adding a model error covariance 
matrix, P'", at the forecast step so that mi; 

p/' = MP"M^-fP"' (27) 

In practice the matrix P'" should be considered as a measure of the variability of the 
noise sequence. This approach has been particularly attractive in the past in view 
of its simplicity and because of the lack of more refined model for the model error 
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dynamics. Note that while is propagated in time and is therefore flow dependent, 
P'" is defined once for all and it is then kept constant. 


3.1 Short Time Extended Kalman Filter - ST-EKF 

We study here the possibility of estimating the model error covariance, P'", on a 
deterministic basis a. The approach uses the formalism on model error dynamics 
outlined in Sect. 2. 

Model error is regarded as a time-correlated process and the short-time evolution 
laws a and (|9]l are used to estimate the bias, b'”, and the model error covariance 
matrix, P'”, respectively. The adoption of the short-time approximation is also legiti¬ 
mated by the sequential nature of the EKF, and an important practical concern is the 
ratio between the duration of the short-time regime and the length of the assimilation 
interval T over which the approximation is used 041 . 

A key issue is the estimation of the two first statistical moments of the tendency 
mismatch, f — f, required in ® and in (|9]l respectively. The problem is addressed as¬ 
suming that a reanalysis dataset of relevant geophysical fields is available and is used 
as a proxy of the nature evolution. Reanalysis programs constitute the best-possible 
estimate of the Earth system over an extended period of time, using an homoge¬ 
neous model and data assimilation procedure, and are of paramount importance in 
climate diagnosis (see e.g. uni). 

Let us suppose to have access to such a reanalysis which includes the analysis, 
x", and the forecast field, so that x{{tj + Tr) = ./^x"(tj), and ly is the assimilation 
interval of the data assimilation scheme used to produce the reanalysis; the suffix r 
stands for reanalysis. Under this assumption the following approximation is made: 

f(x, X) - f(x, = 

x/^(f-f T^)-x"(f) X"(f-f x{{t + T:r)-xf{t + Ty) 5x" 

-—-—-(za) 

'Zf 'Zf 'Zf "Zf 

The difference between the analysis and the forecast, 5x°, is usually referred, in data 
assimilation literature, to as the analysis increment. Erom (l28T l we see that the vector 
of analysis increments can be used to estimate the difference between the model and 
the true tendencies. A similar approach was originally introduced by Leith (1978) 
ca, and it has been used recently to account for model error in data assimilation 

ED. 

Note that the estimate (|28] | neglects the analysis error, so that its accuracy is 
connected to that of the data assimilation algorithm used to produce the reanalysis, 
which is in turn related to the characteristics of the observational network such as 
number, distribution and frequency of the observations. However this error is present 
and acts as an initial condition error, a contribution which is already accounted for 
in the EKE update by the forecast error covariance, As a consequence when (|9]l 
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is used to estimate only the model error component, an overestimation is expected 
that can be overcome by an optimal tuning of the amplitude of b'" and P'”. 

The most straightforward way to estimate the bias due to model error using ( |28] | 
in ©, so that at analysis time it reads: 

b'" = -\/a < 5x“ > - (29) 

Tf 

The bias is then removed from the forecast field before the latter enters the EKF 
analysis update, (l24l i. The scalar term a is a tunable coefficient aimed at optimizing 
the bias size to account for the expected overestimation connected with the use of 
(|28]) . In a similar way the model error contribution to the forecast error covariance 
can be estimated taking the external product of (l28T l after removing the mean and 
reads: ^ 

P'" = a < {5x‘‘- < 5x‘‘ >}{5x«- < K' >}^ > ^ (30) 

V 

We consider now the particular case of parametric error. The forecast error co- 
variance P^, is estimated using the short-time evolution (l2Tli where the correlation 
terms are neglected and the model error covariance, P'" is evolved quadratically 
in the intervals between observations. An additional advantage is that P™ can be 
straightforwardly adapted to different assimilation intervals and for the assimilation 
of asynchronous observations. At analysis times the forecast error bias due to the 
model error, b"\ can be estimated on the basis of the short-time approximation (l20l) : 

b™ =< 5x'” >ft!< > T (31) 

By neglecting the correlation terms and dropping the time dependence for conve¬ 
nience, Eq. (I 2 TI 1 can be rewritten as: 

pf = MP"M^+ < dn„5jxl >z^= MP"M^ + Qt^ = MP"M^ + P"' (32) 

where P" is the analysis error covariance matrix, as estimated at the last analysis 
time, and 

P'" = Qt2 =< SjilJjUo > (33) 

An essential ingredient of the ST-EKE in the case of parametric error is the matrix 
Q: it embeds the information on the model error through the unknown parametric 
error 5X and the parametric functional dependence of the dynamics. In |0 it was 
supposed that some a-priori information on the model error was at disposal and 
could be used to prescribe < ^Mo> and Q then used to compute b'" and P'" required 
by the ST-EKE. The availability of information on the model error, which may come 
in practice from the experience of modelers, is simulated by estimating < 
and Q averaging over a large sample of states on the system’s attractor as. 


b'" = << SjH^ » 

Q = «dnJjHo » 


(34) 

(35) 
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The same assumption is adopted here in the numerical applications with the ST-EKF 
described in Sect. 3.2.1. 

In summary, in the ST-EKF, either in general or in the parameteric error case, 
once b'" and P'” are estimated (with or (ISTT i- dSST l respectively) they are 

then kept constant along the entire assimilation cycle. Model error is thus repeatedly 
corrected in the subspace spanned by the range of P"' where it is supposed to be 
confined. This choice reflects the assumption that the impact of model uncertainty 
on the forecast error does not fluctuate too much along the analysis cycle. Finally, 
in the ST-EKF, the forecast held and error covariance are transformed according to; 




(36) 

P/ = 

^ pf p'« 

(37) 


These new first guess and forecast error covariance, (l3^ and (IJTI i. are then used in 
the EKE analysis formulas (I24li - (l25]) . 


3.1.1 Numerical Results with ST-EKF. Error due to unresolved scale 

We show here numerical results of the ST-EKF for the case of model error aris¬ 
ing from the lack of description of a scale of motion. The case of parametric error 
is numerically tested in Sect. 3.2.1. A standard approach, known in geosciences 
as observation system simulation experiments (OSSE), is adopted here El. This 
experimental setup is based on a twin model configuration in which a trajectory, 
solution of the system taken to represent the actual dynamics, is sampled to produce 
synthetic observations. A second model provides the trajectory that assimilates the 
observations. 

As a prototype of two-scales chaotic dynamics we consider the model introduced 
by 1291, whose equations read: 


dxi , , he ^ r , 

-^ = ixi+\-Xi_2)xi-\-Xi + F -— 36} (38) 

dyji , , he , , 

= -cbyj+i^i(yj+2,i-yi-\,i) - cy;,, -b —x,-, j = {1,..., 10} (39) 

at b 

The model possesses two distinct scales of motion evolving according to (l38l l and 
(EU, respectively. The large/slow scale variable, x,, represents a generic meteoro¬ 
logical variable over a circle at fixed latitude. In both set of equations, the quadratic 
term simulates the advection, the second rhs term the internal dissipation, while 
the constant term in (l38t plays the role of the external forcing. The two scales are 
coupled in such a way that the small/fast scale variables yjj inhibit the larger ones, 
while the opposite occurs for the effect of the variables x, on According to 1291 
the variables yy,- can be taken to represent some convective-scale quantity, while the 
variables x,- favor this convective activity. The model parameters are set as in l29l : 
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c = b = 10, which makes the variables Xj to vary ten times slower than yj i, with 
amplitudes ten times larger, while F = 10 and h = \. With this choice, the dynam¬ 
ics is chaotic. The numerical integration have been performed using a fourth-order 
Runge-Kutta scheme with a time step of 0.0083 units, corresponding to 1 hour of 
simulated time. 

In the experiments the full equations (l38T l - (l39l l. are taken to represent the truth, 
while the model sees only the slow scale and its equations are given by (l3^ with¬ 
out the last term. A network of M = 12 regularly spaced noisy observations of x is 
simulated by sampling the reference true trajectory and adding a Gaussian random 
observation error. We first generate a long record of analysis for the state vector, x, 
which constitutes the reanalysis dataset. The EKF algorithm is run for 10 years with 
assimilation interval t,- = 6 hours, and observation variance set to 5% of the sys¬ 
tem’s climate variance. From this long integration we extract the record of analysis 
increments required in ( |29] | and (l30l l. 

An illustration of the impact of the proposed treatment of the model error is 
given in Fig. [T] which shows a 30 days long assimilation cycle. The upper panel 
displays the true large scale variable xig (blue line), the corresponding estimates 
obtained with the ST-FKF and the FKF without the model error treatment (red and 
yellow lines respectively) and the observations (green marks). The error variance 
of the EKF estimates are shown in the bottom panel. From the top panel we see 
the improvement in the tracking of the true trajectory obtained by implementing the 
proposed model error treatment; this is particularly evident in the proximity of the 
maxima and minima of the true signal. The benefit is further evident by looking at 
the estimated error variance which undergoes a rapid convergence to values close or 
below the observation error. 

A common practical procedure used to account for model error in KF-like and 
ensemble-based schemes, is the multiplicative covariance inflation HI. The forecast 
error covariance matrix is multiplied by a scalar factor and thus inflated while 
keeping its spatial structure unchanged, so that P-^ —(1 + p)P^ before its use in 
the analysis update, (EHi. We have followed the same procedure here and have op¬ 
timized the EKF by tuning the inflation factor p; the results are reported in Fig. 
12 a) which shows the normalized estimation error variance as a function of p. The 
experiments last for 210 days, and the results are averaged in time, after an initial 
transient of 30 days, and over a sample of 100 random initial conditions. The best 
performance is obtained by inflating P^ by 9% of its original amplitude and the es¬ 
timation error variance is about 6% of the system’s climate variance, slightly above 
the observation error variance. Note that when p — 0 Alter divergence occurs in 
some of the 100 experiments. 

We now test the sensitivity of the ST-EKF to the multiplicative coefficient a in 
(|2^ and (l30l l. The results are reported in Fig. |2b), which shows the estimation 
error variance as a function of a. As above the averages are taken over 180 days 
and over the same ensemble of 100 random initial conditions. The important feature 
is the existence of a range of values of a, for which the estimation error is below 
the observation error level. Note that for a = 1, the estimation error is about 4% of 
the climate’s variance, below the observational accuracy. This result highlights the 
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Fig. 1 Top Panel: Model variable xk, for the truth (blue), EKF without model error treatment 
(yellow), EKF with model error treatment (red) and observations (green). Bottom Panel: Estimation 
error variance, normalized with respect to the system’s climate variance, as a function of time. From 

Cl 



Fig. 2 Averaged normalized estimation error variance as a function of (a) the inflation factor p, 
(b) the coefficient a (log scale in the x-axis), and (c) time evolution of the normalized estimation 
error variance for the case p = 0.09 (black) and a = 0.5 (red) (the time running mean is displayed 
with dashed lines). From [7]. 


accuracy of the estimate of P'” despite the simplifying assumptions such as the one 
associated with the correlation between model error and initial condition error and 
the use of the reanalysis field as a proxy of the actual true trajectory. Interestingly, 
the best performance is obtained with a = 0.5, in agreement with the expected 
overestimation connected with the use of (l28l l. 
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In Fig. |3c) we explicitly compare the EKF with the optimal inflation for P^, 
(p = 0.09, P"' = pP^), with the EKF implementing the model error treatment 
through the matrix P'" estimated according to (l30l l and tuned with the optimal values 
of the scalar coefficient a — 0.5. The figure displays the estimation error variance 
as a function of time. Note in particular the ability of the filter, using P'”, to keep 
estimation error small even in correspondence with the two large deviations experi¬ 
enced by the EKF employing the forecast error covariance inflation. 


3.2 Short Time Augmented Extended Kalman Filter - ST-AEKF 

Let us now turn to the state-augmentation approach. In this case we will assume that 
model errors arise from mis-specifications of some parameters, so that the theory 
depicted in Section 2.2 can be used. This view restricts us to parametric errors, but 
it also reflects our limited knowledge of the sub-grid scale processes that are only 
represented through parametrisation schemes for which only a set of parameters is 
accessible. 

A straightforward theory exists for the estimation of the uncertain parameters 
along with the system’s state. The approach, commonly known as state-augmentation 
mi, consists in defining an augmented dynamical system which allocates, along 
with the system’s state, the model parameters to be estimated. The analysis update 
is then applied to this new augmented dynamical system. Our aim here is to use the 
state-augmentation approach in conjunction with the deterministic formulation of 
the model error dynamics. 

The dynamical system (|2^ . the forecast model, is augmented with the model 
parameters, as follows: 


z 


f - 


'xf 


= = 




(40) 


where z = (x, A ) is the augmented state vector. The augmented dynamical system 
includes the dynamical model for the system’s state, and a dynamical model for 
the parameters . In the absence of additional information, a persistence model 
for is usually assumed so that = I and X{^^^ = A“^. Recently, a temporally 
smoothed version of the persistence model has been used in the framework of a 
square root filter Bbl . The state-augmented formulation is also successfully applied 
in the context of the general class of ensemble-based data assimilation procedures 

im. 

By proceeding formally as for Eq. (fT2]) we can write the linearized error evolution 
for the augmented system, in an arbitrary assimilation interval, T = f<;+i — 4, with 
initial condition given by the augmented state analysis error, 5z“ = (5x",5A“) = 

(x"-y,A‘'-A): 
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, rt+T 

5z-/«i(5x^,5A^) = (M5x“ + 2 (41) 

with 5jU“ = (^|;^<i)5A“. The parametric error = 5A^ is constant over the 

assimilation interval in virtue of the assumption = I. Equation (|4TJ describes, 
in the linear approximation, the error evolution in the augmented dynamical sys¬ 
tem drol l. The short-time approximation of the error dynamics (HTI i in the interval T 
reads: 

5zf ^ (M5x" + 5A“) (42) 

As for the standard EKE, by taking the expectation of the product of (HTI) (or 
(|42]) ) with its transpose, we obtain the forecast error covariance matrix, p{, for the 
augmented system: 


P{ =< 5z^5z^^ >= 




(43) 


where the N x N matrix P| is the error covariance of the state estimate x^, P^ is 

f 

the P X P parametric error covariance and P^^ the N x P error correlation matrix 
between the state vector, x, and the vector of parameters A. These correlations are 
essential for the estimation of the parameters. In general one does not have access to 
a direct measurement of the parameters, and information are only obtained through 
observations of the system’s state. As a consequence, at the analysis step, the esti¬ 
mate of the parameters will be updated only if they correlate with the system’s state, 

f 

that is Ki ^ 0. The gain of information coming from the observations is thus spread 
out to the full augmented system phase space. 

Let us define, in analogy with (|4^ . the analysis error covariance matrix for the 
augmented system: 


P 


a 

Z 


na pi 
^ X ^ 

nciT 


xX 


(44) 


where the entries in (l44l i are defined as in ( |4^ but refer now to the analysis step 
after the assimilation of observations. 

By inserting (l42l i into (l43l l. and taking the expectation, we obtain the forecast 
error covariance matrix in the linear and short-time approximation: 


P{ = M < 5x"5x"^ > M^+ < 

[M < 5x"5/x“^ > + < 5/x“5x“^ > M^]t 
= MP“M^ +Q"t2 + [M< 5x“5jti"^ > -b < 5ju"5x"^ > M^]t (45) 
P{ =< 5A“5A“^ > (46) 

p{^ = M < 5x"5A"^ > + < 5 m“5A"^ > T 


(47) 
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Note that ( l45T l is equivalent to (l32l i. except that now the correlations between the 
initial condition and the model error are maintained (last two terms on the r.h.s. 
of (|45]|), and P" and Q“ replace P" and Q. Nevertheless, in contrast to the ST- 
EKF where Q is estimated statistically and then kept fixed, in the ST-AEKF Q" is 
estimated online using the observations. 

The information on the uncertainty in the model parameters is embedded in the 
error covariance P^, a by-product of the assimilation. Using the definition of 
and (|4^ . the matrix Q" can be rewritten as: 


Q" =< >= 


< 








Similarly, the correlation terms in (l45l l can be written according to: 


[M < > M^]t = 


(48) 



Using (|48]) and ( |49] | in (l45l l. the forecast state error covariance vl can be written in 
terms of the state-augmented analysis error covariance matrix at the last observation 
time, according to: 


:MP"M' -f I — |;i,^ 






[MP" 




I A," 


ill 




(50) 


The three terms in ( fSOl l represent the contribution to the forecast state error covari¬ 
ance coming from the analysis error covariance in the system’s state, in the param¬ 
eters and in their correlation respectively. 

By making use of the definition of the model error vector in (l47l i. the forecast 

f 

error correlation matrix P^^ becomes: 


:MP 


xX 


il 


Ia." 




(51) 


Expressions (l46l l. (fSOl l and (fSTT i can be compacted into a single expression: 


P{ = CP“C^ 


(52) 
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with C being the ST-AEKF forward operator defined as: 


C = 



(53) 


where Ip is the P x P identity matrix. 

The short-time bias equation (l20l i is used to estimate the bias in the state forecast, 
x-^, due to parametric error, in analogy with the ST-EKF. This estimate is made 
online using the last innovation of the parameter vector. Assuming furthermore that 
the forecast of the parameter is unbiased, the bias in the state augmented forecast at 
time reads: 

b™ = (54) 

The bias b'” is then removed from the forecast field before the latter is used in the 
analysis update, that is = z-^ — b™, where z^ is the unbiased state augmented 
forecast. 

As for the standard EKE, we need the observation operator linking the model 
to the observed variables. An augmented observation operator is introduced, = 
0] with as in (l2Tt . and its linearization, Hj is now a M x {N + P) matrix in 
which the last P columns contain zeros; the rank deficiency in reflects the lack 
of direct observations of the model parameters. 

The augmented state and covariance update complete the algorithm: 

z" = [I,-K,Hjz^ + K,y“ (55) 

P" = [L-K,HJP{ (56) 

where the vector of observations y“ is the same as in the standard EKE while Ij is 
now the (A + P) x (A + P) identity matrix. The augmented gain matrix is defined 
accordingly: 

K, = P{hJ [H,P{H[ + R] ^ ‘ (57) 

but it is now a (/ + P) xM matrix. 

Equations (l40l i - (l5^ for the forecast step, and (fSSl) - (l57l) for the analysis update 
define the ST-AEKF. The algorithm is closed and self consistent meaning that, once 
it has been initialized, it does not need any external information (such as statistically 
estimated error covariances) and the state, the parameters and the associated error 
covariances are all estimated online using the observations. 

The ST-AEKF is a short-time approximation of the classical augmented EKE, 
the AEKF ifT^ . In essence, the approximation consists of the use of an analytic 
expression for the evolution of the model error component of the forecast error co- 
variance. This evolution law, quadratic for short-time, reflects a generic and intrinsic 
feature of the model error dynamics, connected to the model sensitivity, to perturbed 
parameters and to the degree of dynamical instability. It does not depend on the spe¬ 
cific numerical integration scheme adopted for the evolution of the model state. The 
state error covariance, vl, in the ST-AEKF, is evolved as in the standard AEKF: the 
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propagator M is the product of the individual M,- relative to each time-step within 
the assimilation interval. The difference between the two algorithms is in the time 
propagation of the forecast error covariance associated with the misspecification of 
the parameters, P.A. In the ST-AEKF this is reduced to the evaluation of the off di¬ 
agonal term in the operator C. This term replaces the full linearization of the model 
equations with respect to the estimated parameters, required by the AEKF. In this 
latter case the model equations are linearized with respect to the augmented state, 
(x.A), giving rise to an augmented tangent linear model, Mz. This linearization 
can be particularly involved ll22l . especially in the case of implicit or semi-implicit 
integration schemes such as those often used in NWP applications ETIl . The prop¬ 
agator relative to the entire assimilation interval is then given by the product of the 
individual augmented tangent linear propagator over the single time-steps. As a con¬ 
sequence the cost of evolving the model error covariance in the AEKF grows with 
the assimilation interval. In the ST-AEKF, the use of the short-time approximation 
within the assimilation interval makes straightforward the implementation of the pa¬ 
rameter estimation in the context of a pre-existing EKE, without the need to use an 
augmented tangent linear model during the data assimilation interval. It reduces the 
computational cost with respect to the AEKF, because the propagation of the model 
error component does not depend on the length of the assimilation interval. Never¬ 
theless the simplifications in the setup and the reduction in the computational cost 
are obtained at the price of a decrease in the accuracy with respect to the AEKF. The 
degree of dynamical instabilities along with the length of the assimilation interval, 
are the key factors affecting the accuracy of the ST-AEKF. 


3.2.1 Numerical Results with ST-EKF and ST-AEKF 

Numerical experiments are carried out with two different models. OSSEs are per¬ 
formed first using the Lorenz ’96 ll2^ model used in Sect. 3.1.1, but in its one-scale 
version given by: 

dx' 

-^ = a{xi+i-xi-2)xi-i-I5xi + F, / = {F---,36} (58) 

where the parameter associated with the advection, a, linear dissipation, p and the 
forcing F, are written explicitly. As for the experiments described in Sect. 3.1.1, the 
numerical integration are performed using a fourth-order Runge-Kutta scheme with 
a time step of 0.0083 units, corresponding to 1 hour of simulated time. 

The reference trajectory, representing the true evolution we intend to estimate, 
is given by a solution of (l58l l with parameters = {F‘'', a"', j3"^) = (8,1,1); with 
this choice the model behaves chaotically. A network of M = 18 regularly spaced 
noisy observations is simulated by sampling the reference true trajectory and adding 
a Gaussian random observation error whose variance is set to 5% of the system’s 
climate variance. Model error is simulated by perturbing F, a and j3 with respect 
to their reference true values. Gaussian samples of 100 states and model parameters 
are used to initialize assimilation cycles lasting for 1 year. In all the experiments the 
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initial condition error variance is set to 20% of the system’s climate variance. The 
model parameters are sampled from a Gaussian distribution with mean equal to 
and standard deviation = 25% of 

We compare four configurations of the EKF: (1) standard EKF without model 
error treatment, (2) standard EKE using a perfect model, (3) ST-EKE (Sect. 3.1), 
and (4) ST-AEKE (Sect. 3.2). Recall that in the ST-EKE, model error bias and co- 
variance are estimated according to (l3Tl i and (l33l l with < S^lo> and Q evaluated 
on a statistical basis before the assimilation experiments. The expectation of is 
estimated through; 


<5Mo>=«^Ia(A-A"')» (59) 

and is then used in (IJTT i and (1331) . In ( |59] ) the averages are taken over the same Gaus¬ 
sian sample of initial conditions and parameters used to initialize the data assimila¬ 
tion experiments, using the actual value of the parameter, A^^, as the reference. This 
idealized procedure has been chosen to give the ST-EKE the best-possible statistical 
estimate of the model error in view of its comparison with the more sophisticated 
ST-AEKE. 

Figure |3] shows the analysis error variance as a function of time for the four 
experiments of one year long; the assimilation interval is T = 6 hours. The er¬ 
rors are spatio-temporal average over the ensemble of 100 experiments and over 
the model domain, and normalized with the system’s climate variance. The figure 


T = 6 hrs - = 25% 



Fig. 3 Time averaged analysis error variance as a function of time. Standard EKF without model 
error treatment (black), standard EKF with perfect model (red), ST-EKE (blue) and ST-AEKF 
(green). The error variance is normalized with respect to the system’s climate variance. 


clearly shows the advantage of incorporating a model error treatment: the average 
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error of the ST-EKF is almost half of the corresponding to the standard EKF without 
model error treatment. However using the ST-AEKE the error is further reduced and 
attains a level very close to the perfect model case. 

The benefit of incorporating a parameter estimation procedure in the ST-AEKE 
is displayed in Eig.|4]that shows the time mean analysis error variance for the EKE 
with perfect model and the ST-AEKE (top panel), along with the relative parametric 
errors as a function of time (bottom panel). The time series of the ST-AEKE analysis 
error variance is also superimposed to the time-averages in the top panel. Figure |4] 
reveals that the ST-AEKE is successful in recovering the true parameters. This re¬ 
construction is very effective for the forcing, F, and the advection, a, and at a lesser 
extent for the dissipation, j3. The ability of the ST-AEKE to efficiently exploit the 


Error in the state estimate 



day 


Fig. 4 Top Panel - Time averaged analysis error variance as a function of time: standard EKF with 
perfect model (red) and ST-AEKF (green); time series of the ST-AEKF (black). Bottom Panel - 
Absolute parametric error of the ST-AEKF, relative to the true value . The error variance is 
normalized with respect to the system’s climate variance. 


observations of the system’s state to estimate an uncertain parameter, either multi¬ 
plicative or additive, is evident. Given that the innovation in the parameter, obtained 
via Eq. (fSSl l. is proportional to the cross-covariance forecast error, the accuracy 
of the parameter estimation revealed by Fig. |4] turns out to be an indication of the 
quality of the short-time approximation, (fSTT i. on which the estimate of is based. 

Figure |5] focuses on the comparison between the ST-AEKF and the standard 
AEKF. The experiments are carried out for T = 3, 6 and 12 hours and with 
Ox = 25%X‘''. As above, the results are averaged over an ensemble of 100 exper¬ 
iments, and the observation error variance is 5% of the system’s climate variance. 
The left panels display the quadratic estimation error, while the parametric error is 
given in the panels in the right column; note that the logarithm scale is used in the y- 
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axis. The estimation error relative to the EKF with a perfect model is also displayed 
for reference. 


T = 12hrs 



Fig. 5 Left column: Running mean of the quadratic state estimation error as a function of time; 
EKF with prefect model (red), ST-AEKF (green) and AEKF (blue). Right column: absolute value 
of the parametric error as a function of time for a (red), (blue) and F (green), for ST-AEKF 
(solid lines) and AEKF (dotted lines).From top to bottom T = 12, 6 and 3 hours respectively. The 
errors are averaged over an ensemble of 100 experiments and Oi = 25%A"^. 


We see that as expected the AEKF has a superior skill than the ST-AEKF but for 
T = 3 or 6 hours their performances are very similar. The AEKF shows a marked 
rapidity to reach convergence but the asymptotic error level attained by the two 
filters are practically indistinguishable. On the other hand for T = 12 hours the ST- 
AEKF diverges whereas the AEKF is able to control error growth and maintain the 
estimation error to a low level. We first observe that in all but one cases the para¬ 
metric error in the experiments with the AEKF is lower than for the ST-AEKF, in 
agreement with the observed lower state estimation error. Anyhow when T = 6 or 
3 hours, the asymptotic parametric errors of the two filters are very similar, a re¬ 
markable result considering the approximate evolution law used in the ST-AEKF. 
An important difference is the extreme variability observed in the parametric error 
with the ST-AEKF as compared to the smoothness of the corresponding solutions 
with the AEKF. Note also that when T = 6 hours the ST-AEKF reduces the error in 
the forcing more than the AEKF but the error curves are subject to very large fluc¬ 
tuations. The dissipation, j3, appears as the most difficult parameter to be estimated 
in agreement with what observed in Fig.|4] In summary. Fig. |5] suggests that the ST- 
AEKF may represent a suitable and efficient alternative to the full AEKF when the 
assimilation interval does not exceed the time range of validity of the approximation 
on which the ST-AEKF is based. The results indicate that this limit is between 6 and 
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12 hours given that the ST-AEKF diverges when t = 12 hours. According to the 
theory outlined in ll^ . the short-time regime is related to the inverse of twice the 
largest (in absolute value) Lyapunov exponent of the system. In the Lorenz system 
(fSSl) the largest Lyapunov exponent turns out to be the most negative one, equal to 
0.97 day^*, so that the duration of the short-time regime is estimated to be about 12 
hours, in qualitative agreement with the performance of the ST-AEKF. Finally note 
that the slight deterioration in the filter accuracy is compensated by a reduction in 
both the computational and implementation costs with respect to the AEKF. 

The second model under consideration is an offline version of the operational 
soil model. Interactions between Surface, Biosphere,and Atmosphere (ISBA) ll^ . 
In the experiments that follow, ST-AEKF has been implemented in the presence of 
parametric errors in the Leaf Area Index (LAI) and in the Albedo; more details, 
along with the case of other land surface parameters, can be found in 0. OSSEs 
are performed using the two-layers version of ISBA which describes the evolution 
of soil temperature and moisture contents; the model is available within a surface 
externalized platform (SEDAS; lUO)). The state vector, x = (Tj, 72, W 2 ), contains 

the surface and deep soil temperatures Ts and T 2 and the corresponding water con¬ 
tents Wg and W 2 . The vector X is taken to represent the set of model parameters. A 
detailed description of ISBA can be found in ll^ . 

The forcing data are the same for the truth and the assimilation solutions. They 
consist of 1-hourly air temperature, specific humidity, atmospheric pressure, in¬ 
coming global radiation, incoming long-wave radiation, precipitation rate and wind 
speed relative to the ten summers in the decade 1990-1999 extract from ECMWF 
Re-analysis ERA40 and then dynamically down-scaled to 10 km horizontal reso¬ 
lution over Belgium ifTTll . The fields are then temporally interpolated to get data 
consistent with the time resolution of the integration scheme of ISBA (300 s). In 
this study ISBA is run in one offline single column mode for a 90 day period, and 
the forcing parameters are those relative to the grid point closest to Brussels. An 
one-point soil model has been also used by IIJTI . for parameter estimation using an 
ensemble based assimilation algorithm. 

The simulated observations are 72,,, and RHjm, interpolated between the forcing 
level (Ri20 m) and the surface with the Geleyn’s interpolation scheme ( ifTSl l. at 00, 
06, 12 and 18 UTC. The assimilation interval is T = 6 hours, while the observational 
noise is drawn from a Gaussian, ^(0,R), with zero-mean and covariance given 
by the diagonal matrix R with elements; diag(R) = 

As explained in ED, the observation operatorJ^, relating the state vector to 
the observation includes the model integration. The initial P" and P'” required 
by the EKE are set as diagonal with elements diag{P") = ~ 

(17:2,17:2,10-2, 10-2)Qandflf/flg(P-) ^ ^ (25 x 10-27:2,25 x 

10 - 27 : 2,4 X 10-4,4 X 10-4) 

Parametric errors is introduced by perturbing simultaneously the LAI and the 
albedo. These parameters strongly influence the surface energy balance budget and 


^ The values of CT,,,. and are expressed as soil wetness index SWI — (w — / {wff. — where wj, is the 

volumetric field capacity and Wyt,jif is the wilting point. 
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partitioning, which in turn regulate the circulation patterns and modify the hydro- 
logical processes. For each summer in the period 1990-1999, a reference trajec¬ 
tory is generated by integrating the model with LAI = 1 irP' /n? and albedo = 0.2. 
Around each of these trajectories, Gaussian samples of 100 initial conditions and 
uncertain parameters are used to initialize the assimilation cycles. The initial condi¬ 
tions are sampled from a distribution with standard deviation {ots : , (^wg, <^W 2 ) = 

{5K,5K, 1,1), whereas LAI and the albedo are sampled with standard deviations, 
Cla/ = 0.5 n-p- jm^ and Oaibedo = 0.05 respectively (HU). The initial in the ST- 
AEKF read = 1 {nP jnP'f- and = 10^"*, P" is taken as in the EKF while 

P“ ^ is initially set to zero. 



Fig. 6 RMS e.stimation error in the four state variables for the EKF (red) and the STAEKF (blue). 
The RMS error in the estimate of LAI and Albedo relative to the ST-AEKF are shown in the 
bottom-left/right panels respeetively. From |4l. 


Results are summarized in Eig.|U which shows the RMS Error in the four state 
variable for the EKE and ST-AEKF, along with the RMS Error in LAI and Albedo 
for the ST-AEKF. The progressive parametric error reduction achieved with the ST- 
AEKF is reflected by the systematically lower estimation error in the soil tempera¬ 
ture and water content. At the very initial times, on the order of one week, EKF and 
ST-AEKF have an indistinguishable skill. However, as soon as, the state-parameter 
error correlations in the ST-AEKF augmented forecast error matrix become mature, 
the improvement of the ST-AEKF becomes apparent and it lasts for the entire du¬ 
ration of the experiment. By reducing the parametric error a better guess for the 
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system state can be obtained and this in turn improves the analysis field and again 
the accuracy of the parameter estimate. Moreover, given that this feature is incor¬ 
porated using the short-time formulation 0, the additional computational cost with 
respect to the standard EKF is almost negligible. 


4 Deterministic model error treatment in variational data 
assimilation 

Variational assimilation attempts to solve the smoothing problem of a simultaneous 
fit to all observations distributed within a given interval of interest. We suppose 
therefore that I measurements, ( |2^ . are collected at the discrete times, (fi ,f 2 : 
within a reference time interval T. An priori estimation, x*, of the model initial 
condition is supposed to be available. This is usually referred to as the background 
state, and: 

xo=xi, + Eh (60) 

where e\, represents the background error. 

We search for the trajectory that, on the basis of the background field and accord¬ 
ing to some specified criteria, best fits the observations over the reference period T. 
Besides the observations and the background, the model dynamics itself represents 
a source of additional information to be exploited in the state estimate. The model 
is not perfect, and we assume that an additive error affects the model prediction in 
the form: 

x(f) = .^(xo) 4-5x"’(f) (61) 

Assuming furthermore that all errors are Gaussian and do not correlate with each 
other, the quadratic penalty functional, combining all the information, takes the form 

Gl: 

27= / [ i5x'’'{t')f{P’”)^\,{5x'''{t''))dt'dt'+ 

Jo Jo “ 

f (ek")X'(ek) + ebB-ieb (62) 

A:=l 

The weighting matrices = P(f ,f ),Rk and B have to be regarded as a mea¬ 
sure of our confidence in the model, in the observations and in the background field, 
respectively. In this Gaussian formulation these weights can be chosen to reflect 
the relevant moments of the corresponding Gaussian error distributions. The best- 
fit is defined as the solution, x(f), minimizing the cost-function 7 over the interval 
T. It is known that, under the aforementioned hypothesis of Gaussian errors, x(f) 
corresponds to the maximum likelihood solution and 7 can be used to define a mul¬ 
tivariate distribution of x(f) lfT9l . Note that, in order to minimize 7 all errors have to 
be explicitly written as a function of the trajectory x(f). 
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The variational problem defined by ( l62l i is usually referred to as weak-constraint 
given that the model dynamics is affected by errors iQ). An important particular 
case is the strong-constraint variational assimilation in which the model is assumed 
to be perfect, that is 5x"‘ = 0, ll24l 1^ . In this case the model-error related term 
disappears and the cost-function reads: 

/ 

2Jstrong = L (ek)'^Rk^(ek) + eb B-^eb (63) 

k=l 

The calculus of variations can be used to find the extremum of (l62l i (or (|63]l) and 
leads to the correspondingEuler-Lagrange equations ll24l [^. In the strong-constraint 
case, the requirement that the solution has to follow the dynamics exactly is satisfied 
by appending to (l6^ the model equations as a constraint by using a proper Lagrange 
multiplier field. However the size and complexity of the typical NWP problems is 
such that the Euler-Lagrange equations cannot be practically solved unless drastic 
approximations are introduced. When the dynamics is linear and the amount of ob¬ 
servations is not very large, the Euler-Lagrange equations can be efficiently solved 
with the method of representers 0. An extension of this approach to nonlinear dy¬ 
namics has been proposed in ll44l . Nevertheless, the representers method is far from 
being applicable for realistic high dimensional problems, like the numerical weather 
prediction and an attractive alternative is represented by the descent methods which 
makes use of the gradient vector of the cost-function in an iterative minimization 
procedure ED. This latter approach is used in most of the operational NWP centers 
which employ variational assimilation. Note that in the cost-functions (l62l i model 
error is allowed to be correlated in time, and gives up the double integral in the first 
r.h.s. term. If model error is assumed to be a random uncorrelated noise, only co- 
variances have to be taken into account and the double integral reduces to a single 
integral (to a single summation in the discrete times case). 

The search for the best-fit trajectory by minimizing the associated cost-function 
requires the specification of the weighting matrices. The estimation of the matrices 
P"y, is particularly difficult in realistic NWP applications due to the large size of the 
typical models currently in use. Therefore it becomes crucial to define approaches 
for modeling the matrices P'y, and reduce the number of parameters required for 
their estimation. We will show below how the deterministic, and short-time, model 
error formulation described in Sect. 2.2 can be used to derive P"y, 

We make the conjecture that, as long as the errors in the initial condition 
and in the model parameters are small, the second rhs term of (fT2T i. 5x"'(f) = 
can be used to estimate the model error entering the weak-constraint 
cost-function, and the corresponding correlation matrices P™(f ,f"). In this case, the 
model error dependence on the model state, induces the dependence of model error 
correlation on the correlation time scale of the model variables themselves. By tak¬ 
ing the expectation of the product of the second rhs term of (fT2l i by itself, over an 
ensemble of realizations around a specific trajectory, we obtain an equation for the 
model error correlation matrix: 
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P"'{t\t")=[ dr f t/rM/ , (64) 

Jto Jto ’ ’ ' 

The integral equation ( l64b gives the model error correlation between times t' and t". 
In this form, Eq. ( l64l ) is of little practical use for any realistic non-linear systems. 
A suitable expression can be obtained by considering its short-time approximation 
through a Taylor expansion around (f ,f ) = (foTo)- It can be shown (||5l) that the 
first non-trivial order is quadratic and reads; 

?(/,/')«< 5jUo5/iJ > it - to){t" - to) (65) 

Equation (1651 ) states that the model error correlation between two arbitrary times, t' 
and t , within the short-time regime, is equal to the model error covariance at the 
origin, < S/IgS/XQ >, multiplied by the product of the two time intervals. Naturally 
the accuracy of this approximation is connected on the one hand to the length of 
the reference time period, on the other to the accuracy of the knowledge about the 
error in the parameters needed to estimate < SHo^l^o propose to use the 

short-time law (l65l) as an estimate of the model error correlations in the variational 
assimilation. The resulting algorithm is hereafter referred to as Short-Time-Weak- 
Constraint-4DVar (ST-w4DVar). Besides the fact of being a short-time approxima¬ 
tion, ( l65l l is based on the hypothesis of linear error dynamics. To highlight advan¬ 
tages and drawbacks of its application, we explicitly compare ST-w4DVar with other 
formulations. 


4.1 Numerical Results with ST-w4DVar 

The analysis is carried out in the context of two systems of increasing complexity. 
We first deal with a very simple example of scalar dynamics which is fully inte- 
grable. The variational problem is solved with the technique of representers. The 
simplicity of the dynamics allows us to explicitly solve (l64l i and use it to estimate 
the model error correlations. This ’’full weak-constraint” formulation of the 4DVar 
is evaluated and compared with the ST-w4DVar employing the short-time evolution 
law (l65l l. In addition, a comparison is made with the widely used strong-constraint 
4DVar in which the model is considered as perfect. In the last part of the Section 
we extend the analysis to an idealized nonlinear chaotic system. In this case the 
minimization is made by using an iterative descent method which makes use of 
the cost-function gradient. In this nonlinear context ST-w4DVar is compared to the 
strong-constraint and to a weak-constraint 4DVar in which model error is treated as 
a random uncorrelated noise as it is often assumed in realistic applications. 

Let us consider the simple scalar dynamics; 

A"-; 


r(f) = xoe' 


( 66 ) 
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with X'’' > 0, as our reference. 

Suppose that I noisy observations of the state variable are available at the discrete 
times S [0,T], 1 <k< I: 

yl=^k + ek 

being an additive random noise with variance C7o ( 4 ) = I <k < I, and that a 
background estimate, Xh, of the initial condition, xq, is at our disposal: 


X0=Xh + Eb 


with Eb being the background error with variance (7^. We assume the model is given 
by: 


x(f) =xoe^'. 


We seek for a solution minimizing simultaneously the error associated with all these 
information sources. The quadratic cost function can be written in this case as: 


2J(x) = 


ff 


[xit ) — xoe 


Xt 




(x(f ) —xoe^‘ )dt dt + 


Y^(^o^{yk-^kf + ^b^{^o-Xbf (67) 

A:=l 

The control variable here is the entire trajectory within the assimilation interval T. 
In Eq. ( IFtI ) we have used the fact that the model error bias, 5x™(f), is given by 
x{t) — xqc^^ assuming the model and the control trajectory, x(f), are started from 
the same initial condition xq. Note that xq is itself part of the estimation problem 
through the background term in the cost-function, and that the covariance matrices 
all reduce to scalar, such as p^i^n = p(t ,f ). 

While complete details can be found in 0 we describe here the essential of the 
derivation. The final minimizing solution of (IFTl i is found using the technique of 
representer and reads: 


x{t)=xbe^'+'^^krk{t)=x^{t)+'^hrk{t) 0<t<T (68) 

k=l k=l 


The I functions, rk{t), are the representers given by: 

rk{t) = rk{0)e^‘+ [ p^At)ak{t')dt 

Jo " 


\<k<I 


(69) 


subject to r^(0) = Jq ak{t)e^At, \<k<I, while the adjoint representers satisfy: 

ak{t) = 8{t — tk) ^<k<I (70) 

subject to ak{T) = 0, 1 < k < /. The coefficients, j3^, are given by: 

j3 = (S + (J„%)-'d 


( 71 ) 
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with d the innovation vector, d = (/j — xf, ), S the / x / matrix {S)ij = 

ri{tj), and 1^/ the / x / identity matrix. The coefficients are then inserted in dfiSl l to 
obtain the final solution. 

In the derivation of the general solution dfiST l (with the coefficients (ITTl i). we have 
not specified the model error correlations p^{t ,t ); the particular choice adopted 
characterizes the formulations we aim to compare. Our first choice consists in eval¬ 
uating the model error correlations through (l64l i. By inserting = ■^SX, with 
/(x) = Ax, and the fundamental matrix, associated with the dy¬ 

namics dfifil) . we get: 


p^{t /) =< (xodXf > \'t" (72) 

where the expectation, <>, is an average over a sample of initial conditions and 
parametric errors. Expression (|7^ can now be inserted into (I69]l- (l70l i. to obtain the 
I representer functions in this case: 

r^(f)=e^('+'‘)[<(xo5A)2>4f-f(J^2] \<k<M (73) 

The representers dT^ are then inserted into ( ItTI ) to obtain the coefficients for the so¬ 
lution, x(f), which is finally obtained through dfiSl l. This solution is hereafter referred 
to as the full weak-constraint. 

The same derivation is now repeated with the model error weights given by the 
short-time approximation (l65l l. By substituting 5/i = into (l65l l. we obtain: 

p^(t ,t ) =< (xo5A)^ > 11 (74) 

Once (|74| is inserted into ( l69b - iTOl) the representer solutions become: 

rk{t) = < (xo5A)2 >tkt \<k<I (75) 

The representer functions are then introduced into ItTI i and (l68l) to obtain the so¬ 
lution, x(f), during the reference period T. The solution based on (ITST l is the ST- 
w4DVar. 

The strong-constraint solution is derived by invoking the continuity of the so¬ 
lution (|7^ . or (iTSl l. with respect to the model error weights. The strong-constraint 
solution is obtained in the limit 5A —0, and reads: 

rk{t) = l<k<I (76) 

The three solutions based respectively on (|7^ . (fTSl) and (l76l) are compared in 
Fig. [T] Simulated noisy observations sampled from a Gaussian distribution around 
a solution of (l66l) . are distributed every 5 time units over an assimilation interval 
r = 50 time units. Different regimes of motion are considered by varying the true 
parameter A"”. The results displayed in Fig. |2]are averages over 10^ initial condi¬ 
tions and parametric model errors, around xq = 2 and A"” respectively. The initial 
conditions are sampled from a Gaussian distribution with standard deviation <7;, = 1, 
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while the model parameter. A, is sampled by a Gaussian distribution with standard 
deviation \ AX \ = \X'’' — A |; the observation error standard deviation is Go = 0.5. 

Figure |7] shows the mean quadratic estimation error, as a function of time, during 
the assimilation period T. The different panels refer to experiments with different 
parameter for the truth 0.01 < A"- < 0 .03, while the parametric error relative to the 
true value is set to AXjX'’’ — 50%. The three lines refer to the full weak-constraint 
(dashed line), ST-w4DVar (continuous line) and the strong-constraint (dotted line) 
solutions respectively. The bottom right panel summarizes the results and shows the 
mean error, averaged also in time, as a function of A"" for the weak-constraint so¬ 
lutions only. As expected the full weak-constraint solution performs systematically 
better than any other approach. ST-w4DVar successfully outperforms the strong- 
constraint case, particularly at the beginning and end of the assimilation interval. 
The last plot displays the increase of total error of this solution as a function of A^'^. 
To understand this dependence, one must recall that the duration of the short-time 
regime in a chaotic system is bounded by the inverse of the largest amplitude Lya¬ 
punov exponent ll^ . For the scalar unstable case considered here, this role is played 
by the parameter A"”. The increase of the total error of the short-time approximated 
weak-constraint as a function of V'' reflects the progressive decrease of the accuracy 
of the short-time approximation for this fixed data assimilation interval, T. 

The accuracy of the ST-w4DVar in relation to the level of instability of the dy¬ 
namics, is further summarized in Fig. [8] where the difference between the mean 
quadratic error of this solution and the full weak-constraint one, is plotted as a func¬ 
tion of the adimensional parameter rA"", with 10 < T <60 and 0.0100 < A"" < 
0.0275. In all the experiments AA/A"" = 50%. Remarkably all curves are superim¬ 
posed, a clear indication that the accuracy of the analysis depends essentially on the 
product of the instability of the system and the data assimilation interval. 

We turn now to the case of a nonlinear dynamics. We adopt here the widely used 
Lorenz 3-variable convective system Il28l , whose equations read; 

dx 

dy 

= px-y-xz ill) 

dt 

dz „ 

with A = ((7,p,j3) = (10,28, |). OSSEs are performed with a solution of (l77l) rep¬ 
resenting the reference dynamics from which observations are sampled. The estima¬ 
tion is based on observations of the entire system’s state (i.e. observation operator 
equal to the identity 3x3 matrix), distributed within a given assimilation inter¬ 
val and affected by an uncorrelated Gaussian error with covariance R. The model 
dynamics is given by (l77l i with a modified set of parameters. The numerical inte¬ 
grations are carried out with a second order Runge-Kutta scheme with a time-step 
equal to 0.01 adimensional time units. 

The variational cost function can be written, according to (l62] |. as: 
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l' = .0100 l' = .0125 




Fig. 7 Mean quadratic estimation error as a function of time for variational assimilation with sys¬ 
tem f66t . The panels refer to experiments with different A"'. The bottom-right panel shows the 
mean quadratic eiror, for the weak-constraint solutions only, averaged also over the assimilation 
interval T as a function of A"'. Strong-constraint solution (dotted line), full weak-constraint solu¬ 
tion (dashed line), short-time approximated weak-constraint solution (continuous line). From (^. 


L L 

27(xo,xi,...,xz,) = 

i=i;=i 




l)) + 


^(y^-Jf'(xA.))^R '(y^-Jf^(xi-))-F(xfc-xo)^B '(x^-xq) (78) 

k=\ 

We have assumed the assimilation interval T has been discretized over L time steps 
of equal length. At. 

The control variable for the minimization is the series of the model state x, at 
each time-step in the interval T. The minimizing solution is obtained by using a 
descent iterative method which makes use of the cost-function gradient with respect 
to x/, 0 <i <L. This latter reads: 
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Fig. 8 Difference between the mean quadratic error of the short-time approximated and the full 
weak-constraint solution, for system 1^, as a function of for assimilation period 10 < T’ < 

60 and for different values of In all the experiments = 50%. From (5]. 


V,„y = -HjR-'(yg-jr(xo)) 

-Mji (P'f y)-' (xj- - (x;-i))] - B-\xh - xo) i = 0 

y=i 

V,,J = -HfR-i (y“ - jr(x,)) - M^,+i [£ (x,- - 

7=1 

L 

+ L ^7j i^j - 1)) 1 < / < L - 1 (79) 

7=1 

= -H^R-i(y“-jr(x^))+f (P;^,,.)^‘(x7-^(x,-i)) i = L 

7=1 

The gradient ( l79b is derived assuming that observations are available at each 
time step f,, 0 < t < L. In the usual case of sparse observations the term proportional 
to the innovation disappears from the gradient with respect to the state vector at a 








Deterministic treatment of model error in geophysical data assimilation 


33 


time when observations are not present. Note furthermore that, if the model error is 
treated as an uncorrelated noise, the corresponding term in the cost-function reduces 
to a single summation over the time-steps weighted by the inverse of the model error 
covariances. The cost-function gradient modifies accordingly and the summation 
over all time-steps disappears ll42l . 

The cost-function (fTsT l and its gradient (|79]) define the discrete weak-constraint 
variational problem. The ST-w4DVar consists in (fTsT l and (1?^ with the model er¬ 
ror correlations P”f estimated using the short-time approximation (IhST i that, in this 
discrete case, reads; 

p;.",. =< 5Ho5^ll > ijAt^ (80) 

The invariant term < >, which is here a 3 x 3 symmetric matrix, is as¬ 

sumed known a-priori and estimated by accumulating statistics on the model attrac¬ 
tor, so that < 5Hq5hI >=« 5Hq5hI » and perturbing randomly each of the 
three parameters <7, p and j3, with respect to the canonical values and with a stan¬ 
dard deviation \AX\. The ST-w4DVar is compared with the weak-contraint 4DVar 
with uncorrelated model error formulation and with the strong-constraint 4DVar; in 
this latter case the model error term disappears from the cost-function ( l78b and the 
gradient is computed with respect to the initial condition only BTIl . 

The assumption of uncorrelated model error is done often in applications. It is 
particularly attractive in view of the consequent reduction of the computational cost 
associated with the minimization procedure. Model error covariance are commonly 
modeled as proportional to the background matrix, so that P'" = aB. Figure|9]shows 
the mean quadratic error as a function of the tuning parameter. Results are averaged 
over an ensemble of 50 initial conditions and parametric model error; the obser¬ 
vation and assimilation interval are set to 2 and 8 time-steps respectively. Strong 
constraint 4DVar (green) and ST-w4DVar (red) do not depend on a and are there¬ 
fore horizontal lines in the panel. Weak constraint 4DVar with uncorrelated model 
error (blue) shows, as expected, a marked dependence on the model error covari¬ 
ance amplitude. The blue line with squared marks refers to an experiment where the 
model error is treated as an uncorrelated noise but the spatial covariances at observ¬ 
ing times are estimated using the short-time approximation. Comparing this curve 
with the blue line relative to weak constraint 4DVar with uncorrelated model error 
and P'" = aB allow for evaluating the impact of neglecting the time correlation and 
of using an incorrect spatial covariance. 

The uncorrelated noise formulation (solid line with no marks) never reaches the 
accuracy of the ST-w4DVar. Note furthermore that for small a it almost reaches 
the same error level as the strong-constraint 4DVar where the model is assumed to 
be perfect. By further increasing a over a = 10^ (not shown) the error reaches a 
plateau whose value is controlled by the observation error level. When the spatial 
covariances are estimated as in the short-time weak-constraint the performance is 
generally improved, although for large a, the estimate P'" = aB gives very close 
skill and the improvement in correspondence with the best-possible a is only mi¬ 
nor. This suggests that the degradation of the uncorrelated noise formulation over 
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AX/X"" = 1 0% 



Fig. 9 Mean quadratic estimation error as a function of the tuning parameter a multiplying the 
model error covariance in the weak-constraint 4DVar with the uncorrelated noise assumption (see 
text for details). The dynamics is given by system Cj}. ST-w4DVar (red), strong-constraint 4DVar 
(green), uncorrelated noise weak-constraint 4DVar (blue) and uncorrelated noise weak-constraint 
4DVar with spatial covariance as in the short-time approximated weak-constraint (blue with red 
marks). From 


the short-time weak-constraint is mainly the consequence of neglecting the time 
correlation and only to a small extent to the use of an incorrect spatial covariance. 


5 Discussion 

Data assimilation schemes are usually assuming the uncorrelated nature of model 
uncertainties. This choice is indeed legitimate as a first approximation when initial 
condition errors dominate model errors. Due to the large increase of measurement 
data availability and quality, this view should be reassessed. However one promi¬ 
nent difficulty in dealing with model errors is the wide variety of potential sources 
of uncertainties, going from parametric errors up to the absence of description of 
some dynamical processes. But recently a stream of works, caEilEa, provided 
important insights into the dynamics of deterministic model uncertainties, and from 
which generic mechanisms of growth were disentangled. In particular, it was shown 
that the mean square error associated with the presence of deterministic model errors 
is growing quadratically in time at short lead times. These insights now open new 
avenues in the description of data assimilation schemes, as it has been demonstrated 
in the present chapter. 

First, the deterministic approach was applied in the context of the Extended 
Kalman Filter, for both the classical state estimation scheme and its augmented ver- 
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sion (state and parameters). It has been demonstrated that these new schemes are in¬ 
deed most valuable when dealing with deterministic model error sources. They pro¬ 
vide a large improvement in the state estimation (in particular with the ST-AEKF) 
not only in the context of idealized settings (Lorenz’ system) but as also for realistic 
applications as shown by the results obtained with the offline version of an oper¬ 
ational soil model (ISBA). Second the same idea was investigated in the context 
of four dimensional variational assimilation for which a weak-constrained frame¬ 
work was adopted. In this case model error cross-correlations were considered as 
quantities depending quadratically on time, implying a time dependent weighting 
of the model error terms during the assimilation period. This approach also led to 
important improvements as compared to more drastic assumptions like the time in¬ 
dependence of model error source terms, or the absence of such sources. 

The approaches proposed in this chapter rely on an important assumption, the 
deterministic nature of model uncertainties. As alluded in Section 2, the proposed 
general setting could also be extended to independent random noises, provided the 
appropriate temporal variations of the bias (Eq. 5) and covariances (Eq. 7) are used 
(see the remark at the end of Section 2). Still the covariance will be dependent on 
time (linearly) (see e.g. EsJ) and will for instance affect differently the weights of 
the weak constrained four dimensional variational assimilation. Besides these tech¬ 
nical aspects, it remains difficult to know what the exact nature of the sources of 
model errors is. A realistic view would be that the fast time scale processes - like 
turbulence in the surface boundary layer - could be considered as random compo¬ 
nents, while parameterization errors in the larger scale description of the stability of 
the atmosphere is a deterministic process. This leads us to consider model errors at a 
certain scale of description as a deterministic component plus a random component, 
and the user is left with the delicate question of evaluating whether the dominant 
sources present in his problem are of one kind or the other. 

This work has mostly tackled this problem in the context of simple idealized 
dynamical systems. These encouraging results need however to be confronted with 
more operational problems. In line with the results obtained with soil model, this 
problem is currently under investigation in the context of the use of the ST-AEKF 
for an online version of ISBA coupled with ALARO at the Royal Meteorological 
Institute of Belgium. A relevant methodological development will be the incorpo¬ 
ration of the deterministic model error treatment in the context of the ensemble 
based schemes, as illustrated in or ll38l . Finally, it is worth mentioning that the 
deterministic model error dynamics has been recently used in a new drift correc¬ 
tion procedure in the context of interannual-to-decadal predictions 0. These recent 
applications illustrate the usefulness of the deterministic approach and should be 
further extended to a wider range of applications in which model error is present. 
In particular in the context of coupled atmosphere-ocean systems where multiple 
scales of motion are present and model error often originates at the level of the 
coupling. 
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